library(tidyverse)
library(cowplot)
library(patchwork)
library(daymetr)
library(raster)
library(spatstat)
library(KrigR)
library(lubridate)
library(sf)
library(viridis)
library(terra)
library(XPolaris)
library(MetBrewer)Soil data
Downloading and wragling soil data
Packages
About XPolaris
XPolaris has been removed from CRAN. To install the GitHub available version:
if (!require(XPolaris)) {
if (!require(devtools)) {
install.packages("devtools")
}
devtools::install_github("lhmrosso/XPolaris")
}White mold data
Loading white mold data. We are going to use the field coordinates (latitude and longitude) to extract the data from the rasters of soil variables.
If you have not downloaded the data necessary to run these analysis, please refer to getting started section before runnig the code below
wm_data = read.csv("data_white-mold/WhiteMoldSurveyWrangledData.csv")Removing missing coordinates
There are some missing coordinates in the dataset. Here we remove them.
wm_data2 = wm_data %>%
filter(!is.na(latitude))Soil variables
Locations
Setting up the coordinated of each quadrat we need to download the soil data
# exkansas
ny_locations = data.frame(ID = c("NY1","NY2","NY3","NY4","NY5","NY6","NY7","NY8"),
lat = c(42, 42, 42, 42, 43, 43, 43, 43),
long = c(-77, -78,-79, -80,-77, -78.0,-79.0, -80))
xplot(locations = ny_locations)+
geom_point(data = wm_data2 , size = 0.3,
aes(longitude,latitude))+
coord_map( xlim = c(-81,-75),
ylim = c(41.5,44.5))Coordinate system already present. Adding new coordinate system, which will
replace the existing one.
# max(wm_data2$longitude)Download soil data
ny_images <- ximages(locations = ny_locations,
statistics = c('mean'),
variables = c('ph','om','clay',"sand","silt","bd", "hb","n","alpha","ksat","lambda","theta_r","theta_s"),
layersdepths = c('0_5'),
localPath = file.path("soil_images"))Here we read the metadata of the downloaded data. It containg the information regarding the directory of each downloaded file.
Merge images
Organic matter
We are going to filter only the lines contain information regarding Organic matter
om_df_images = ny_images %>%
filter(variables == "om")As you can notice in the next figure, the raster objects are downloaded in separate raster of 01 degree size
# read all files for organic matter
om_stack = lapply(om_df_images$local_file, brick)
par(mfrow = c(2,4))
for(i in 1:length(om_stack)){
plot(om_stack[[i]])
}Therefore, we use the function merge() to create a single raster object covering the whole study region
om_ny_raster =merge(om_stack[[1]],
om_stack[[2]],
om_stack[[3]],
om_stack[[4]],
om_stack[[5]],
om_stack[[6]],
om_stack[[7]],
om_stack[[8]])
plot(exp(om_ny_raster))Automatization
Instead of doing the above step for each variable by hand, here we created a function to do all the steps automatically.
get_soil_var = function(var, data){
data_filtered = data %>%
filter(variables == var)
stack_list = lapply(data_filtered$local_file, brick)
merged_raster = merge(stack_list[[1]],
stack_list[[2]],
stack_list[[3]],
stack_list[[4]],
stack_list[[5]],
stack_list[[6]],
stack_list[[7]],
stack_list[[8]]
)
return(merged_raster)
}
soil_vars = c('ph','om','clay',"sand","silt","bd", "hb","n","alpha","ksat","lambda","theta_r","theta_s")
selected_vars = c('ph','om','clay',"sand","silt","bd","theta_r","theta_s")soil_variables_list = lapply(soil_vars,get_soil_var, data = ny_images)
names(soil_variables_list) = soil_vars
saveRDS(soil_variables_list, "soil_images/list_soil_variables_raster.rds")aggre_var_list = lapply(soil_variables_list, aggregate,fact=30)
saveRDS(aggre_var_list, "soil_images/list_soil_variables_raster_aggregated.rds")Plot soil maps
NY shape file
Loading New York state shape file.
ny_shape1 = read_sf("shape_files/cugir-007865/cugir-007865/cty036.shp")Cropping raster files
We use the function lappy() to crop all variables’ rasters using the NY shape file as a mask.
aggre_var_list2 = lapply(aggre_var_list, mask, ny_shape1)Then here we create a function for ploting the soil maps.
actual_var_names = c("Soil pH in water", "Soil organic matter","Clay","Sand","Silt","Bulk density","Residual soil water content","Saturated soil water content")
actual_var_symbol = c("pH", "OM","Clay","Sand","Silt","BD","\u03B8r","\u03B8s")
actual_var_units = c("", "(%)","(%)","(%)","(%)","(g/cm³)","(m³/m³)","(m³/m³)")
plot_gg_raster = function(X,raster, var){
# actual_var_names[X]
if(var[X] == "om"){xx=1}else{xx = 0}
as.data.frame(raster[[var[X]]], xy = T) %>%
filter(layer !="NaN", x< -76.8) %>%
mutate(layer = case_when(xx ==1~ exp(layer),
xx ==0~ layer)) %>%
ggplot(aes())+
geom_raster(aes(x, y, fill = layer))+
scale_fill_viridis(option ="B",guide = guide_colorbar(barwidth = 0.2, barheight =5 ))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
linewidth =0.2,
alpha = 0.5,
color = "white")+
# geom_point(data = wm_data2 , size = 0.1,color = "white",
# aes(longitude,latitude))+
# coord_quickmap(xlim = c(-80,-76.8), ylim = c(42,43.35))+
coord_sf(xlim = c(-80,-76.8), ylim = c(42,43.35))+
theme_map()+
labs(title =paste(" ",actual_var_names[X]),
fill = paste(actual_var_symbol[X],actual_var_units[X]))
}
# selected_vars
# aggre_var_list[[1]]
# plot_gg_raster(1,aggre_var_list2, var = selected_vars[1] )as.data.frame(aggre_var_list2$ph, xy = T) %>%
filter(layer !="NaN", x< -76.8) %>%
# mutate(layer = case_when(xx ==1~ exp(layer),
# xx ==0~ layer)) %>%
ggplot(aes())+
geom_raster(aes(x, y, fill = layer))+
scale_fill_viridis(option ="B",guide = guide_colorbar(barwidth = 0.2, barheight =5 ))+
geom_sf(data = ny_shape1,
# aes(x=long, y = lat, group = group),
fill= NA,
linewidth =0.2,
alpha = 0.5,
color = "white")Combo soil maps
Here we plot all maps into a single combo fiure
do.call(patchwork::wrap_plots, lapply(X =1:length(selected_vars) , FUN =plot_gg_raster, raster = aggre_var_list2, var = selected_vars))+
plot_layout(ncol = 2)+
plot_annotation(tag_levels = "A")&
theme(legend.position = "right",
legend.text = element_text(size = 5),
legend.title = element_text(size = 5),
plot.title = element_text(size = 7, face = "bold"))ggsave("figs/maps/soil_maps.png", dpi = 900, height = 7, width = 7, bg = "white")Extract variables to location
wm_data2_uni = wm_data2 %>%
group_by(subject) %>%
slice(1L)Selecting coordinate columns from the white mold data set
coords<-data.frame(lon=wm_data2_uni$longitude, lat=wm_data2_uni$latitude)
coordinates(coords)<-c("lon","lat")Extracting variable from the original merged raster (30 meters resolution)
df1 = lapply(soil_variables_list, extract, coords@coords)
as.data.frame(df1) %>%
mutate(subject =wm_data2_uni$subject) %>%
cbind(longitude=wm_data2_uni$longitude,
latitude=wm_data2_uni$latitude) %>%
write.csv("soil_images/extracted_soil_data.csv",row.names = F)Session Info
sessionInfo()R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=Portuguese_Brazil.utf8 LC_CTYPE=Portuguese_Brazil.utf8
[3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C
[5] LC_TIME=Portuguese_Brazil.utf8
time zone: America/Sao_Paulo
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] knitr_1.48 MetBrewer_0.2.0 XPolaris_1.0.2
[4] terra_1.8-5 viridis_0.6.5 viridisLite_0.4.2
[7] sf_1.0-19 KrigR_0.9.4 spatstat_3.3-0
[10] spatstat.linnet_3.2-3 spatstat.model_3.3-3 rpart_4.1.23
[13] spatstat.explore_3.3-3 nlme_3.1-164 spatstat.random_3.3-2
[16] spatstat.geom_3.3-4 spatstat.univar_3.1-1 spatstat.data_3.1-4
[19] raster_3.6-30 sp_2.1-4 daymetr_1.7.1
[22] patchwork_1.3.0 cowplot_1.1.3 lubridate_1.9.3
[25] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[28] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
[31] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] DBI_1.2.3 pbapply_1.7-2 deldir_2.0-4
[4] gridExtra_2.3 rlang_1.1.4 magrittr_2.0.3
[7] e1071_1.7-16 compiler_4.4.1 mgcv_1.9-1
[10] systemfonts_1.1.0 maps_3.4.2.1 vctrs_0.6.5
[13] ecmwfr_2.0.3 crayon_1.5.3 pkgconfig_2.0.3
[16] fastmap_1.2.0 labeling_0.4.3 utf8_1.2.4
[19] rmarkdown_2.28 tzdb_0.4.0 ragg_1.3.3
[22] automap_1.1-12 xfun_0.48 cachem_1.1.0
[25] jsonlite_1.8.9 progress_1.2.3 goftest_1.2-3
[28] reshape_0.8.9 spatstat.utils_3.1-1 prettyunits_1.2.0
[31] parallel_4.4.1 R6_2.5.1 stringi_1.8.4
[34] stars_0.6-7 Rcpp_1.0.13 iterators_1.0.14
[37] tensor_1.5 zoo_1.8-12 snow_0.4-4
[40] FNN_1.1.4.1 Matrix_1.7-0 splines_4.4.1
[43] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.17.0
[46] abind_1.4-8 yaml_2.3.10 codetools_0.2-20
[49] curl_6.0.1 lattice_0.22-6 intervals_0.15.5
[52] plyr_1.8.9 withr_3.0.2 evaluate_1.0.1
[55] units_0.8-5 proxy_0.4-27 polyclip_1.10-7
[58] xts_0.14.1 pillar_1.9.0 KernSmooth_2.23-24
[61] renv_1.1.2 foreach_1.5.2 ncdf4_1.23
[64] generics_0.1.3 spacetime_1.3-2 hms_1.1.3
[67] munsell_0.5.1 scales_1.3.0 class_7.3-22
[70] glue_1.8.0 mapproj_1.2.11 tools_4.4.1
[73] grid_4.4.1 colorspace_2.1-1 cli_3.6.3
[76] spatstat.sparse_3.1-0 gstat_2.1-2 textshaping_0.4.0
[79] fansi_1.0.6 doSNOW_1.0.20 gtable_0.3.5
[82] digest_0.6.37 classInt_0.4-10 htmlwidgets_1.6.4
[85] farver_2.1.2 memoise_2.0.1 htmltools_0.5.8.1
[88] lifecycle_1.0.4 httr_1.4.7